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We investigate the finite temperature momentum distribution of bosonic noncondensate particles 
inside a 3D optical lattice near the superfiuid to Mott insulator transition point, treating the quan- 
tum fiuctuation and thermal fluctuation effects on equal footing. We explicitly address the different 
momentum (q) dependence of quasi-particle distribution resulted from thermal and quantum ori- 
gins; the former scales as |q|~^ and hence is dominant in the small momentum region, while the 
later scales as |ql~^ and hence dominant in the large momentum limit. Analytic and semi-analytic 
results are derived, providing a unique method to determine various properties inside the optical 
lattice, including temperature, condensate density, coherent length and/or single particle gap etc. 
Our results also agree with the scaling theory of a quantum XY model near the transition point. 
Experimental implication of the TOP measurement is also discussed. 
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I. INTRODUCTION: 

The experimental realization of superfluid (SF) to 
Mott insulator (MI) transition of ultracold atoms 
in an optical lattice has open a new area of strongly cor- 
relatedphysics and lead to many applications to other 
fields It is generally believed that when the size of 
the expanding atom cloud is much larger than the initial 
size (i.e. a long-time fiight after switching off the trap- 
ping potential) and if the interaction effect during the 
expansion can be neglected [3| , the TOF image can be in- 
terpreted as a momentum distribution function of the ini- 
tial atom cloud inside the optical lattice. As a result, the 
sharp interference peaks at zero momentum and Bragg 
momentum can be understood as a signature of a Bose- 
Einstein condensation (or superfluidity) while the wider 
and smaller hump around the peak is then attributed 
to the non-condensate particles inside the optical lat- 
tice. Different from the superfluidity measurement of 
condensate via vortices in a single parabolic confinement 
potential (5|, the time-of- fiight (TOF) absorption image 
(i.e. the " bi-modal" structure with a " sharp" interference 
peak in the TOF image) is so far the only experimen- 
tal measurement to determine the superfluid/condensate 
density inside the optical lattice [11^. 

However, behind the scenario of "bi-modal" structure 
of the interference peak, there is an important assump- 
tion: only two length scales (or momentum scales) are 
relevant in the long wavelength limit and their differ- 
ence can be distinguishable Q : one is associate with the 
size of condensate and the other is associate with the 
thermal wavelength of non-condensate particles. In fact, 
from experiment point of view, only after the later part 
(non-condensate particles) has been identified and sub- 
tracted from the TOF image spectrum, one can study 
the properties (i.e. the sharpness, width, or condensate 
fraction) of the former (condensate) part, which is the 
key player to determine the many-body phase diagram 
and has been extensively discussed in recent theoretical 
and experimental groups jlili^Iililililll ■ To the 



best of our knowledge, so far all the experimental data 
of these non-condensate particles were fitted by a Gaus- 
sian type distribution function at finite momentum and 
the obtained fitting parameter is interpreted as the tem- 
perature of bosons inside the optical lattice 0. This 
approach may be justified for weakly interacting bosons 
at finite temperature, but it cannot be reliable near the 
SF-MI transition point, where the quantum depletion is 
known significantly enhanced. In other words, near the 
quantum transition point, there will be at least three rel- 
evant length scales in the interference peak: The first 
one is the condensate size as mentioned above, the sec- 
ond one is the thermal wavelength, and the third one is 
the healing length, which is associate with the interac- 
tion effect or quantum depletion. (The lattice constant 
and the size of Wannier function can be assumed not 
relevant when investigating a single peak of the interfer- 
ence pattern.) To the best of our knowledge, there is no 
useful analytical form of the momentum distribution to 
distinguish the two different contributions (thermal and 
quantum) of the non-condensate particles near the SF- 
MI transition point. Exact numerical simulation cannot 
distinguish these two contribution either (9llT3| . Without 
a justified theory to describe the momentum distribution 
of non-condensate particles, the measurement and/or de- 
termination of the condensate fraction in the interference 
peak also becomes questionable. Solving such problem 
and providing a useful theoretical framework is the mo- 
tivation and the theme of this work. 

In this paper, we apply the three-state effective model 
developed by Altman et al. to calculate the finite 
temperature momentum distribution of quasi-particles 
near the SF-MI transition, where the quantum fiuctu- 
ation can be comparable or even stronger than the fi- 
nite temperature effect. Our results show that (1) in the 
superfluid regime, the momentum distribution of ther- 
mal excited quasi-particles is different from the one of 
quantum depleted quasi-particles: the former diverges as 
|q|~^ in the small momentum (q) limit, while the lat- 
ter diverges as |q|~^- In other words, the noncondensate 
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particle is dominated by thermal particles in small mo- 
mentum regime while it is by quantum depletion in the 
larger momentum regime. (2) In the weak interacting 
regime and in the Mott insulator regime, we obtained an 
analytical result for the momentum distribution at finite 
temperature, providing the theoretical fitting equations 
for the TOF image. The obtained equations can be used 
to calculate various physical properties, like the conden- 
sate fraction, temperature, and single particle gap etc. 
inside the optical lattice. Note that here we just concen- 
trate on the momentum distribution of particles inside 
the optical lattice, and assume that this distribution can 
be observed in a long TOF experiment when the inter- 
action effects during the expansion is negligible For 
simplicity, we will just consider the uniform system with- 
out the inhomogeneous trapping potential. 

This paper is organized as following: In Section HIl we 
briefly review the system Hamiltonian and the connec- 
tion between momentum distribution inside the optical 
lattice and the momentum distribution in free space. The 
later can be understood as a result of TOF absorption 
image if assuming no interaction effect during expansion 
and neglecting the flnite size effect [l^. In Section Hill we 
briefly introduce the three state model and its meanfield 
phase diagram obtained near the SF-MI transition point, 
both at zero temperature and finite temperature regime. 
In Section [TVl we first study the momentum distribution 
calculated in the weakly interacting regime, where one 
can apply Bogoliubov theory to study the general prop- 
erties of thermal and quantum depleted quasi-particles. 
In Section IVIl we applied the meanfield solution of the 
three state model to the Mott insulator phase and obtain 
some analytic form of momentum distribution function 
at finite temperature regime. Results near the quantum 
critical point (QCP) of the SF-MI transition are inves- 
tigated in Section IVIII In Section |Vl we further apply 
the three-state model to study the momentum distribu- 
tion in the superfluid side near QCP, both below and 
above Tc. We then study and compare the momentum 
distribution near SF-MI transition by using the scaling 
theory of quantum XY model in Section [Villi and then 
discuss some issues related to experimental observation. 
We finally summarize our results in Section [X] All the 
supplementary materials about the details of calculation 
are in the Appendices. 



II. SYSTEM HAMILTONIAN AND 
MOMENTUM DISTRIBUTION 



where or is the bosonic field operator at site R, and 
Tifi — a^aj^ is the density operator. J is the single par- 
ticle tunneling amplitude. U is onsite interaction energy 
and fi is the chemical potential. For simplicity, we ne- 
glect the harmonic trapping potential throughout this 
paper. Within this single band approximation, we can 
express a bosonic field operator (^'(r)) in lattice momen- 
tum expansion (see Appendix]^: 'l'(r) — X^k '^k'I'k(r), 
where Ok = J2ii'^'R.^~^^^ with the lattice vol- 

ume, fl = L'^ (we assume L lattice sites in each direc- 
tion). ^k(r) is the Bloch wavefunction function and Y^'^ 
is the summation of lattice momentum over the first Bril- 
liouin Zone (BZ). As a result, the free space momentum 
distribution, N(q), can be calculated by mapping the 
field operator into the plane wave basis (see Appendix 

m, i.e. A^(q) = AEL"(k)|«'k(q)f: where ^^{q) is 
the Fourier transform of Bloch wavefunction 5'k(r), and 

n(k) = i'^u.O'k) is the momentum distribution of the lat- 
tice momentum, k. 

When the lattice potential is sufficiently strong, the 
lowest band Wannier function can be well-approximated 
by a Gaussian function with width, a. As a result, 
the Bloch wavefunction can be easily calculated to be 
(see Appendix El): |^'k(q)P = /(k - q)|w(q)P, where 



/(q) ^ 



-|q| 



_ (270; 



and 



sin(L<}n/2) 
sin(g„/2) 

l^^(q)l'^ = ^;S72-e~''*' " ■ The momentum distribution 
measured in the TOF experiment is then just a col- 
umn integration of N(q) along the z direction, leading to 
(see Appendix Rl for details) Af_L(q_L) = J^^dqzN{q) cx 

2 

n Ski fc. ''^(k) ^k_L(*lJ-) upto an overall constant. 
Note that to derive above expression, we have used the 
fact that az\kz\ <C 1 in most regime of the first BZ. 
Finally, we note that when the system size is much 
larger than the lattice constant, /(q) becomes a sharp 
peaked function at each reciprocal lattice, Gj_, with 
width ^ 2tt/L. As a result, the TOF image, iVj^(qj^), 
can be further simplified by approximating /(q_L) to be 
a delta function at Gj^, and become (see Appendix \K\i 



N 



'^(q^) (X i^'n(q^,A:,)|w;(q^)p^<5(q^-G^) 



Gj 

CX n_L(q_L)|w(q_L)P as q_L 0, 



(2) 



where ?T.j_(q_L) = X]fc ''^(qj-j k^) is the momentum dis- 
tribution integrated over the z direction. 



In this paper, we consider bosonic atoms loaded in a 
3D square optical lattice with lattice constant, a. When 
the lattice strength is sufficiently large, the system can 
be well described by a single band Hubbard model: 



H = -J 



(R,R') 



u 



— ^nR(nR - 1) -/i^nR,(l) 



R 



R 



III. THREE-STATE MODEL AND FINITE 
TEMPERATURE PHASE DIAGRAM 

The phase diagram of the single band Hubbard model 
has been extensively studied in the literature in the 
last ten years. In 2D and 3D systems with a square 
lattice, meanfield phase diagram has been shown to 
be a very good approximation for the phase boundary 
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FIG. 1: Meanfield phase diagrams of Bose Hubbard model in 
a 3D square lattice for no — 1: (a) finite temperature phase 
diagram and (b) quantum phase diagram. The MI lob in (b) 
is plotted according to meanfield approximation, with the lob 
tip bending down from U /2. The horizontal dashed line in (b) 
indicates the particle-hole symmetric line. The finite temper- 
ature phase diagram in (a) is calculated from Eq. (|16|) . and 
the critical point at T = is right the same as the tip posi- 
tion of the MI lob, as indicated by the vertical dash-dot line. 
Filled dots are the calculated Tc and the thin line is plotted to 
guide the eyes. The thin dashed line above Tc is the line of the 
constant peak width, ko = 0.0087r/a, which is the half-width 
of a Lorenzian fitting function for the momentum distribu- 
tion. The thick dashed lines is the single particle gap energy 
in the Mott insulator phase. ■ ■ • A*'5, are different regimes 
of normal state. SFl and SF2 are superfluid states in weakly 
interacting and strongly interacting regime respectively. 



TlillllH. The meanfield study of the SF-MI phase 



transition can be summarized by the two phase diagrams 
in Fig. [TJa) and (b), where the quantum critical point, 
Uc, is defined at the tip of the MI lob (with average no 
particle per site) in (b). The finite temperature phase 
diagram in (a) is obtained by fixing the fiUing fraction 
to be an integer, (tir) = hq — 1. In order to study 
the momentum distribution in the single band Hubbard 
model, we separate the study into two regimes: one is 
the weakly interacting regime of SF phase, and the other 
is the strongly interacting regime, where superfluid order 
parameter is very small or even zero when entering the 
MI regime. For the convenience of later discussion, below 
we briefly outline the underlying meanfleld theories for 
these two different regimes of phase diagram. 

In the weakly interacting regime, i.e. away from the 
quantum phase transition point in the region 7V1 and 
SFl of Fig. [Ha), we could apply Bogoliubov-Hartree- 




FIG. 2: Probabilities of different number occupation, n, of 
each lattice site in the superfluid state of unit filling (no — 
(c'ri^r) = 1) at T = (see the text). The three solid lines are 
for n = 1, n = 0, and n = 2 respectively from top to bottom. 
The dashed line is for n — 3, and results for larger number 
occupation is not visible in this parameter regime. 



Fock-Popov approximation to study the quasi-particle 
excitation and quantum depletion effects both for T > Tc 
and T < Tc- Since this approach has been very well 
known in the literature [ITJ, we will just outline the re- 
sults and implication to the momentum distribution in 
the next section. However, in the second (strongly in- 
teracting) regime, we adopt the truncated phase space 
method developed by Altman et al. (Ref. [1^) to study 
the superfluid phase near the critical point and the Mott 
insulator phase. Such truncated phase method is justifled 
by using the fact that in the strongly interacting regime, 
number fluctuation per lattice site can be so small that 
only three states with number of particles uq, uq + I, and 
no — 1, are relevant in the low energy regime. In order to 
verify the validity of such truncated phase space approx- 
imation in the superfluid state, we use Gutzwiller trial 
wavefunction, l^'triai) = Hr {J2'^=o fn\n)-R), to calculate 
the probability of higher number occupation at zero tem- 
perature, where |n)R is the wavefunction for n particles 
at site R and /„ is the amplitude of the wavefunction of 
n particles at each lattice site. By minimizing the total 
energy of Eq. ([T]) with a unit flUing fraction (no = 1) 
in average, in Fig. [21 we flnd the probability to have 
more than two particles per site is very small (~ 1%, 
the dashed line) even for U/J ~ 22, away from the SF- 
MI transition point, Uc/ J = 34.97. Since the transition 
temperature is also of the same order as the tunneling 
strength, J, it is reasonable to expect that such trun- 
cated phase space method should be justified in a wide 
range near the SF-MI transition point and at finite tem- 
perature not too far away from Tc. In this paper, we will 
name this truncated phase space scheme to be a three- 
state model. Since details of the theory have been derived 
in Ref. |T^, below we just listed some important results 
for completeness and for the convenience of later discus- 
sion. Details of the momentum distribution calculation 
are shown in the Appendix [B] and [Cl 

In the three-state model, one first truncates the phase 
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space to no — 1, no and uq + 1 particles per site (with 
no being the average fiUing fraction), by defining three 
bosonic operators, ^, to be j^|0) = |no + a)^ — 

(at^)"0^ 



= |0) for a 



0, ±1. Within this truncated 



V(no+Q)! ' 

space, a bosonic creation operator can be expressed to 



be a 



then rewrite the original Hamiltonian of Bose-Hubbard 
model in terms of such new field operators. The com- 
pleteness of such truncated phase space is ensured by 
applying a constrain: X]a=-i ^1 R^a r ~ ^- order to 
apply a variational meanfield study (i.e. Gutzwiller-type 
ground state), we make a unitary transform to a coherent 
state basis: tji — U ■ 6r, where Tr = [ij, ^, t\ tl-^ j^]-^ 

and 5r = [^o r, r, ^2 r]"'"- transform matrix, 

U, depends on some variational parameters (see Ref. 
[l6| and Appendix [B|, which are determined by min- 
imizing the variational energy, (^mf|^^|^a/f)j with 
I^Pmp) = nR^oRl^^) being the trial wavefunction at 
T = 0. Such meanfield calculation can give a mean- 
field quantum phase boundary between SF phase and 
MI phase, as shown in Fig. [Ijb). The tip position in (b) 

/no)^, and 



can be found to be at Uc = zJ{y/no + 1 
Mc = ('^0 — ^)U—^ (here z = 6 is the number of nearest 
neighboring sites in a 3D square lattice). Fixing the aver- 
age density to be an integer, no, the chemical potential at 
T = changes as a function of U as the dashed line of Fig. 
[Ijb). Defining Sfi = fi— (no — as the chemical poten- 
tial deviation from the middle of MI lob at ?7 = 00, such 
particle-hole symmetric (i.e. integer filling) line can be 
obtained to be: Sfi = —j [zJ + U{y/no + 1 + y/no)~^] in 
the SF side and — —zJ/2 in the MI side. The finite 
temperature phase diagram in (a) is plotted along such 
integer filling line. 



To calculate the quasi-particle excitation and the as- 
sociated quantum depletion energy, Altman et al. (l6j 
further used Holstein-Primakoff bosons to eliminate 60, R 
and 6j in order to obtain an effective Hamiltonian to 
the quadratic order of quasi-particles, 61/2, r- In the SF 
phase, the obtained effective Hamiltonian can be further 
diagonalized by using a canonical transformation in mo- 
mentum space (i.e. a generalized Bogoliubov transfor- 
mation), b]^ ^ M ■ Pi^, where 5k = [^i,k' ^2,k' ^i,k' 4,k]^ 

and /3k = [/5s,k'/3m,k'/3l,k'/3m,k]'^ (see Appendix U. 
Here A4 is the transformation matrix, and Ps/mM ^re 
the quasi-particle field operators for the sound mode 
and the massive mode respectively. In the long wave- 
length limit, their dispersion becomes es(k) = Ws|k| and 
e„i(k) — Am + 'k'^ /2m* , where Vg is sound velocity, is 
the mass gap and m* is the effective mass of the massive 
excitation. Using ^ — {Uc — U)/J to measure the dis- 
tance from the critical point, Uc, we find that for integer 



filling (along the dashed line of Fig. [Ub)): 

_ TT^/no{l + no)C/z 



(3) 



^no(l + ?io)(l + 27io) - 2no(l + np) ^ ^ ^^^^^ 



(4) 



near the critical point (i.e. ^ ^ 1). Similar to the ordi- 
nary Bogoliubov transformation in the single component 
case, the quantum depleted particles can be also derived 
from such canonical transformation, which is very signif- 
icant when near the quantum critical point, Uc- We note 
that recent calculation 18] on the spectral weights of 
these excitations show that the massive excitation in the 
superfluid regime has a rather short life time compared 
to the low energy phonon excitations. This result, how- 
ever, cannot be correctly captured within the three-state 
model used in our recent paper. 

In the Mott state, the ground state can be described 
by the three-state model even better, due to the strong 
reduction of number fluctuations. Within the meanfield 
approximation, the MI ground state is still described by 
a product state with integer number of particle per site, 

I'^mf) = IIr^UIO) = IIr^o.rIO)- a result, the 
obtained effective Hamiltonian can be easily diagonalized 
by another canonical transformation [l6j : 



(5) 



where Pp/h^k is the field operators for particle/hole exci- 
tation. Here A{k) = cosh(i:)k/2) and B{k) = sinh(Dk/2) 
are coefficients with 



tanh(£'k) 



-2eo(k)Y/»-o(no 1) 
i7-eo(k)(2no + l) ' 



(6) 



and eo(k) = 2 J ^^^-^ cos(fcQa). The particle and hole 
excitation dispersion are given by 

ep(h)(k) = T[eo(k)/2-H<5A^]-fc:;(k), (7) 

where (5/i is the chemical potential measured from (no — 
i)J7, eo(k) 2JX;a=i cos(A:a), and 

c^(k) = iv/C^2_f/,^(k)(4„^ + 2) + eo(k)2. (8) 

At zero temperature, the particle and hole excitations 
have gaps: 



A = ep(0) + £,,(0) = yJU"^ - zt/J(4no + 2) + (zJ)2. (9) 



Near the critical point, we have A 



2 (no (no 



\)fl^yfzmj ^0{U'^^'^) with bU = U-Uc. In the strong 
interacting limit, we have A = U + 2zJ yJnp{nQ + 1) -|- 
0{U~^). In finite temperature, the chemical potential 
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((5/i(T)) is given by fixing the number of particles to be 
no, i.e. number of particles are the same as number of 
holes: ^ J^w [/^(ep.k) - /^(eh.k)] = 0. More details can 
be found in Appendix [Cl 

Using the meanfield ground states and quasi-particle 
excitations described above, in the rest of this paper we 
will study the momentum distribution at finite temper- 
ature. Our results should be able to give a quantitative 
prediction for the study of SF-MI quantum phase transi- 
tion in the experimentally relevant regime. 



IV. MOMENTUM DISTRIBUTION IN THE 
WEAKLY INTERACTING SUPERFLUID PHASE 

We first study the momentum distribution in the 
weakly interacting regime by using the Bogoliubov- 
Hartree-Fock- Popov approximation [13. We could sep- 
arate the calculation in two different regimes of temper- 
ature: (i) T > Tc, (2) T < Tc. They correspond to 
regions, Nl, and SFl respectively in Fig. [T[a). In the 
normal state regime, it has been shown jl7| that the 
self-energy shift to the single particle dispersion is just 
a constant, ek = Cb^k + 2noJ7, where eb,k = 2zJ — eo(k) is 
the single particle band energy measured from the bot- 
tom (k = 0) and Uq is the average filling fraction. As 
a result, the chemical potential is also shifted by a con- 
stant and the transition temperature, T^, is the same 
as the noninteracting system, giving by (set ks = 1) 
"0 = hT.'ki^^'"'''^^^'^'' - 1)"^- For a 3D square lattice, 
we have TJJ = 5.86 for uq = 1 and Tc/J = 10.29 for 
no — 2. When T is very close to Tc from above, the mo- 
mentum distribution can be well approximated to be a 

in the long 



40r 



Lorentzian function: n'^^(k) 



wavelength limit, where k^i^a 



(|k|2+fc^^Q2 



2nQU\/J 



1.073(5T/Tc)i/2 is ti^g thermal wavevector and 5T = T- 
Tc- When T Tc, the thermal momentum kf^ ""/a, 
and therefore we could calculate A^_L(q_L) by integrat- 
ing kz from — CX3 to CX3 as shown in Eq. We ob- 
tain 7VNi(q_L) = _=^^^==|u;(q_L)|2, where w{qj_) oc 

I |2 2 

g-iq^l is a broad Gaussian function due to the small 
onsite Wannier function of width a < a within the sin- 
gle band approximation. As pointed out in Ref. j^, the 
TOF image obtained from N^-^ (q_L ) can therefore have a 
"sharp peak" even when T > Tc- As a result, such TOF 
image can be easily mis-interpreted as the characteris- 
tic feature of condensation. This simple analysis shows 
the importance for an quantitative analysis for the TOF 
image in the optical lattice. 

When the temperature is below Tc (i.e. region SFl 
of Fig. [ija)), we can treat Oq = oq as a c-number (i.e. 
condensate at k = 0) and apply the Bogoliubov-Hartree- 
Fock-Popov approximation to diagonalizc the quadratic 
order of the resulting effective Hamiltonian. Following 
the standard approach [l7| . we can obtain the phonon 
dispersion: eB(k) 
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FIG. 3: Momentum distribution of non-condensate particles 
for (a) T/J = 0.2 and (b) T/J = 0.02 in region SFl of Fig. 
[TJa). Thin solid and dashed lines are nqn^(k) and nyf^(k) 
respectively, and thick solid lines are their sum. Here we use 
U/ J — 0.5, no = 1 and ky — kz = 0. Approximated results of 
Eq. (llOf) cannot be distinguished with the total distribution 
(thick line). 



momentum distribution can be divided into two parts: 



,SF1 



(k) = n^Qjj^(5(k)-t-n^Jf^(k), where nfon is the number 



density of condensate particles and n^J (k) 



-.SFl 

'■th 



(k)- 



nqjj^(k) is the momentum distribution of non-condensate 



particles. Here nf^^ (k) 



.,SF1 



(k) 



_ £b(k)+raoC/o 



on 

esW) — /s(es(k)) and 
i are the contribution from 



et,(k)+rtoC/o 



^/e6(k)(eb(k) + 2noUo), and the 



"qn V"/ 2eB(k) 2 

thermal excitation and quantum depletion respectively. 
/b(x) = (e^/-^ — 1)~^ is Bose-Einstein distribution func- 
tion. Since the condensate part is a structure-less delta 
function in momentum space, here we concentrate on the 
momentum distribution of non-condensate particles. In 
Fig. [3l we show typical results of 7i^^^(k) and ^^^^^(k) 
for comparison. We note that, even without adding the 
Wannier function, n^jf^(k) and nqjj^(k) can have differ- 
ent shape in both long and short wavelength limit: as 
shown in Fig. [3jb), nn^"^(k) is dominated by thermal 
excitation in the long wavelength limit, while dominated 
by quantum depletion in the large momentum regime. 
This feature also makes it easy to be mis-interpreted as 
a "bimodal structure" of condensate, although we do not 
include the condensate particles yet. This situation will 
become more serious and crucial when near the quantum 
critical point, where the number of condensate particles 
can be small while the "bimodal structure" still exists 
due to the different momentum distribution of thermal 
and quantum excitation. 

To further study the underlying physics in this regime, 
we note that there are two momentum scales involved in 
above expression: one is thermal momentum, kf^^a = 
\/T I J and the other is quantum momentum, k^^a = 

V'^i;^t/72^. There are three regimes we need to iden- 
tify: (i) In weak interaction limit, we have fcq^^ ^ 
k^^ and hence the most relevant momentum regime is 
^qn"'^ |k| ^ ^fif^j dominated by the thermal excita- 
tion. As a result, using eB(k) '-^ e6(k) for |k| > k^^ , 

we find 7i^P(k) ~ n^^'^^{^^) ^ ■^]'kp~ ^^^^ region, (ii) 
In extremely low temperature region, we have fcfjf ^ ^ 
/Cq^^, and the most relevant momentum regime becomes 
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kf^^ |k| <^ fcqn^: dominated by the quantum deple- 
tion instead. Neglecting the thermal contribution in this 

region, we find n^^^(k) ~ nq^^(k) ~ (iii) ^^'^ i^" 

termediate regime (U ~ T), we have fc^jf ^ ~ k^^-^ and the 
most relevant momentum regime is from |k| < fc^jf ^, fcqn ^• 
The the leading contribution of non-condensate particles 
can be obtained by expanding in the long wavelength 
limit: 



,SF1 



(k) 



2|k|2a2 



Co - C2|kpa^ 



Co 



2|k|2a2 l + (c2/co)|k|2a2 



(10) 



where 



Co 



24 



+ 



1 



(fcSFl)2 6(fcfFl)2 2 



and 



C2 



1 (fcfif^)^ 

32 (fcSFl)4-90 (fcSFl)6 



1 (4^')' (MTl! 

480 



12(fcfFia)2- 



(12) 



Our numerical calculation shows that Eq. (jlOp is a very 
good approximation for the full Bogoliubov result (< 1%) 
for T/ J > 0.1 and noU/J < 10, covering most superfluid 
regime that current experiment can access. After column 
integration over kz (we can extend the integration range 
to (— oo, oo) if only A:^if/q„ ^ tt), we find 



7Vfl(q^) 



(k- 



SFl^i2 
th ) 



c2/2 



4|q-L| v/c2(co + C2|qj.|2) 




in this regime. Eq. (fT3|) can be used to fit the experi- 
mental data to obtain the temperature and average filling 
fraction inside the optical lattice (We note that the over- 
all amplitude is arbitrary and therefore only two param- 
eters in Eq. (|13p are needed). In above analysis, we did 
not discuss the structure of condensate part, which is sup- 
posed to be a sharp peak at q = with width tt/L. The 
details of the condensate shape may be strongly affected 
by the interaction effect in the beginning of expansion. 
On the other hand, our analysis on the non-condensate 
particles above provides a quantitative and analytic equa- 
tion for experimentalists to study the thermal/quantum 
tails for q_L ^ 0. These parts is much less affected by 
interaction during the expansion and therefore can be a 
much more reliable method to extract the temperature 
inside optical lattice. 



V. MOMENTUM DISTRIBUTION IN THE 
STRONGLY INTERACTING SUPERFLUID 
PHASE NEAR SF-MI TRANSITION POINT 

When the interaction U is large enough, the quantum 
depletion of the condensate particles becomes significant. 



and therefore the Bogoliubov approach cannot be justi- 
fied, especially when close to the SF-MI transition point. 
More precisely, for the case of unit filling, uq = 1, the Bo- 
goliubov approach used above is justified when U <ti J 
(i.e. fc|f = y/noU/2J < 1), while the SF-MI transi- 
tion occurs at U ^ J. In the later regime, the number 
fiuctuation at each lattice size is strongly reduced and the 
phase fiuctuation becomes enhanced. In this regime, we 
can apply the three-state model to study the SF phase 
as briefly described in Section Hill To calculate the mo- 
mentum distribution in the optical lattice, (al^ci-i^), we 
need to apply a series of transformation to diagonalize 
the quantum fluctuation on top of the meanfield result 
of the three-state model. Details of the calculation is 
shown in Appendix|B] Here we just shown the final result 



of the noncondensate part: n^c (k) 



SF2fi^\ _ ^SF2 



(k)- 



tSF2 
'■qn 



(k), 



where the thermal (n^j^ (k)) and the quantum depleted 
(riq,'j2(k)) parts are respectively (see Eq. (jBlSp ): 



<^(k) 



(5n(k)2 + 533(k)2)/(e,^k) 
+ (522(k)2-F544(k)2)/(e™ 

533(k)2 544(k)^ 



(14) 
(15) 



where Sij (k) is the matrix element of a 4 x 4 matrix 
defined in Eq. (|Bl2|) . The TOF image, ^j_(qj.), can 
be also calculated directly from Eq. ^ based on above 
result. Unfortunately, the analytic expression for the mo- 
mentum distribution is too complicated to be expressed 
even in the long wavelength limit. Therefore in this sec- 
tion we will just show the numerical results in various 
parameter regime. 

Before showing the numerical result of momentum dis- 
tribution, it is instructive to mention that above result 
can be used to calculate the interaction effect of the SF 
transition temperature, Tc(t/), (as shown in Fig. [Ha)), 
by requiring the conservation of total number of particles: 



SF2 



(k) = no. 



(16) 



In other words, when T < Tc, the difference between the 
total number of particles and the non-condensate parti- 
cles can be interpreted as the condensate particles, pro- 
viding the superfluidity of the system. When the in- 
teraction is stronger and/or the temperature is larger, 
the condensate density becomes smaller so that the non- 
condensate particles, contributed both from the thermal 
excitation and the quantum depletion, becomes domi- 
nant. Although we just consider the quantum fluctuation 
to the quadratic order within the three-state model, the 
obtained transition temperature (shown in Fig. Hl^a)) 
is qualitatively consistent with the general picture, de- 
creasing dramatically near the SF-MI transition point at 
T = 0. 

When the temperature is above Tc, the three-state 
model is still justified if only the temperature is much 
smaller than the on-site interaction U (i.e. small number 
fluctuation). Although the disappearance of the order 
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FIG. 4: Momentum distribution of non-condensate particles 
in 3D optical lattice with U = 34J and no = 1 for different 
temperature: T/J = 0.1, 1, 1.5, and 1.7 respectively for (a)- 
(d). The solid and dashed lines are quantum deplation and 
thermal particles, and the thick solid lines are the total. Here 
Tc = 1.2J and the upward arrows in (a) and (b) indicate the 
presence of condensate particles at k = 0. 



VI. MOMENTUM DISTRIBUTION IN MOTT 
INSULATOR PHASE 

Now we move to the study of finite temperature mo- 
mentum distribution in the Mott insulator regime (in re- 
gions N4, and N5 of Fig. [ija)), where the particle number 
fluctuation at each lattice site is strongly reduced. As 
mentioned in Section IIIIl in the three-state model [l^ 
the quasi-particle excitation can be calculated analyti- 
cally, leading to a powerful tool to study the momen- 
tum distribution analytically. After some straightfor- 
ward but tedius calculation (details are in Appendix [Cjl . 
the full momentum distribution can be also contributed 
from thermal part and quantum deplation part as before: 
nf^{\i) = <i(k) + n^^i(k), where 

<i(k) = -5n + /B(ep.k)(V^^j;TTA(k)- V^S(k))' 
+/s(e/>,k) (\^^S(k) - V^A(k))' , (17) 
n^„i(k) = (V^I^i3(k)- VH^A(k))' (18) 



parameter (or condensation) makes the meanfield ap- 
proach of the three-state model not well-justified, but 
quadratic expansion of the elementary excitaion (i.e. 
6i/2.r) in the effective Hmailtonian can be still a good ap- 
proximation if the number fluctuation is still small due 
to the strong interaction in the low temperature reime 
(i.e. U ^ T > Tc). As a result, we can still calculate the 
momentum distribution above Tc by tuning the chemi- 
cal potential away from the value inside the superfluid 
phase, in order to conserve the total number of particles. 
In fact, such deviation of the chemical potential above 
Tc automatically leads to the disappearance of the Gold- 
stein mode and then renormalizes the effective mass of 
the quasi-particle excitation. 

In Fig. |3]we show the numerical results of the momen- 
tum distribution in the optical lattice for the thermal 
and quantum deplated particles in this regime. Differ- 
ent from the results of weakly interacting regime (SFl), 
the quantum deplation can still dominate the momentum 
distribution of non-condensate particle near the critical 
point, Uc, at low temperature regime (Fig. H^a)). When 
the temperature is approaching Tc from below, however, 
the thermal excited particles become dominant in the mo- 
mentum distribution (Fig. IDJb)). Similar to the situation 
of weakly interacting regime discuss above, momentum 
distribution of both quantum deplated and thermal ex- 
cited particles are divergent in the long wavelength limit 
for T < Tc, making the distinguishment between conden- 
sate particles and non-condensate particles highly non- 
trivial, at least in a uniform system. When the temper- 
ature is above Tc, the contribution of condensate parti- 
cles disappears and the momentum distribution becomes 
a Gaussian-like function with the width increased as a 
function of temperature. 



as shown in Eq. (jC24p . Here A(k) and -B(k) have been 
defined in Eqs. ([5]) and ([6]) above. Sn a constant density 
shift due to thermal excitation, and can be calculated 
from Eq. (jC7P (see more details in Appendix [Cj. 

Since the single particle tunneling, J, is much smaller 
than onsite interaction, U, in this regime, there are only 
two temperature limit we need to consider: (i) in the 
low temperature regime (U ^ T, J, i.e. region N4 of Fig. 
[Ija)), and (ii) in the high teperature regime {T ^ U ^ J, 
i.e. region N5 of Fig. [l^a)). From the calculation 
shown in Appendix [Cl in the low temperature region 
{U ^ T,J), the chemical potential can be calculated 



to be (see Eq. (|C15|) : 



^o(^) 



, where 



/o(a;) is the modified Bessel function. As a result, we find 
the momentum distribution of quantum deplated parti- 
cles and the thermal excited particles can be calculated 
explicitly to the leading order of J /U and T /U (see Eq. 
(lck|) ): 



<(k) 



no 



2?T,o(l -I- no] 
U 



eo(k). 



(19) 



-rh'(k) 



5ie 



(no + l).o(k)/T ^ ^^gno.o(k)/T _ (20) 



where 



(1 



92 = J^oe "'""e are small constants in 



3/2 



5n — (2no -I- 



this low temperature limit. 

density shift due to thermal excitation as shown in Eq. 
(IC7|) . The TOF image, A^_L(qj_), can be also calculated 
analytically in the leading order of large U to be (see 



is a constant 
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FIG. 5: Momentum distribution of non-condensate particles 
in the MI regime, (a) and (b) are for T / J = 1, and (c) and 
(d) are for T/J = 10. The left pannels are disribution along 
kx axis with ky = kz — 0, while the right pannels are distri- 
bution along the diagonal direction, i.e. k = fc(l, 1, 1). Thin 
dashed and thin solid lines are the thermal and quantum con- 
tribution respectively, and the thick solid lines are the total 
distribution. Here we use U/ J = 40, and other parameters 
are the same as Fig. O Note that the thermal contribution 
in (d) becomes negative in the large momentum regime due 
to the constant density shift, &n. The total distribution, of 
course, is always positive definite. 



Eqs. (IC251) and (IC26l) : 



4no(l +no)J, . . ... 
no H — (cos(qa;) -h cos(qy)) 



where = ^^/^(M!^) and = g2h(' t >■ 

One interesting thing in the above close form ex- 
pression is that the momentum distribution is different 
from the pure meanfield result (A^j^(qj^) = no) by ad- 
ditional contribution from both thermal fluctuation of 
quasi-particles and the quantum deplation. The ther- 
mal excited particles can contribute a very narrow peak 
when T < J, because e"o'o(«-'°)/'^ oc e2"f''™>'('?-«)/^ ~ 
g2noJ/T(i-g^a /2) ^ showing a characteristic momentum 
scale, •y/T/noJf?, in the long wavelength limit. Although 
such narrow peak structure on top of a uniform momen- 
tum distribution looks surprising, but it is consistent with 
the sharpe peak structure of a normal state just above 
the superfluid Tc in the weakly interacting regime Q . 

Now we consider the high temperature regime when 
T ~ ;7 > J, as the region N5 in Fig. [TJ In this 
regime, the chemical potential correction, (5/.t, can be 
shown to be very small, to the second order of J/T: 
^ = -3(2«o+i) coth(:g.) (^)^ as explicitly calculated 
in Eq. (jC22[) . This result indicated that the chemical 



k(q±)l'(2i) 



' 2,1 ng ' 



potential, ^, is almost the same as the value at the mid- 
dle of MI lob in the large U limit even in the system 
of low filling fraction (no ^ 1), reflecting the fact that 
the particle and hole excitations are almost the same 
(particle- hole symmetry) in this high temperature limit. 
As a result, to the leading order of small single particle 
tunneling, J, the momentum distribution can be calcu- 
lated to be (see Eq. (IC27P and Appendix [Cl for details): 



,N5 



(k) 



,N5 
Hh 



(k) 



,N5 



(k), where 



<'(k) 
<^(k) 



no -I- 2no(l -|- no) 



fo(k) 

U 



= hi 



eo(k) eo(k) 



4no(l+no) 

qU/2T_1 1 



(22) 
(23) 

and 



Here we define the two constants, hi 

/i2 — - — (^^u/2T_iy • Note that the contribution of 

Sn is cancelled in the regime, see Eq. (jC27|) . It is easy 
to see that in this regime, the thermal excitation con- 
tributes the same momentum distribution (oc eo(k)) as 
the one contributed by quantum deplation. This is very 
different from the result in the low temperature limit (see 



-,N4 
hh 



(k) above). This is because, in the low temperature 
limit, the thermal wavelength can be very long, leading 
to a small momentum distribution of the thermal ex- 
cited quasi-particles. However, when the temperature is 
increased, the thermal wavelength becomes shorter and 
even the same order as the lattice constant, and hence 
the momentum distribution is broad in momentum space. 
The TOF image obtained by integrating over the com- 
ponent can be also calculated easily (see Eqs. (|C28|) and 
(EMI): 



^±(q±) = 



no -I- Hi 



eo(q_L) 



k(q±)P (24) 



where Hi 

(l+2no+2nl)e"'^'^ 
(eU/^Tliy 



2na{l+no),J 
U 



eU/2T_i 



VII. MOMENTUM DISTRIBUTION AT 
QUANTUM CRITICAL POINT 

When near the quatum phase transition point (say 
N3 of Fig. [1]), the meanfield treatment becomes un- 
justified due to large fluctuations. However, we find 
that it is still possible to extract some useful informa- 
tion about the momentum distribution in the large no 
limit, where the meanfield approximation becomes jusi- 
fied again, since the transition becomes more classical 
like. We can start from the analytic results in MI state by 
taking U Uc- Using the fact that A(k) = cosh(Dk/2) 
and i3(k) — sinh(Z?k/2) with being defined in Eq. 
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(|6|), we find that in the large Uq limit: 



(k) = feiep^w) {V^k^A{k) - ^B{k)y 



6n 



+ (1 + fBi^hM)) (xAI^S(k) - V^Aik)) 

noil + fsiepM) + /s(e/.,k))(^(k) - B{k)f. 

(25) 



We also have 



{A{k)-B{k)f = 



lU- £o(k)(2no + 1 - 2^no{no + 1)) 
U - eo(k)(2no + 1 + 2y/MmTTj) 



1 7k, 

u 



(26) 



where u = U/inozJ and 7k = z^^J2s^^^'^ with S be- 
ing the direction to the nearest neightboring sites [l6l |. 
In other words, for 3D square lattice, eo(k) = J7k = 
2J^^ cos(fcQ,a). Within the same approximation (no S> 

1), we have ep^k ^ e^.k = (C^/2) yl — Tk/", ^^^^ therefore 
the momentum distribution becomes 



nr (l + 2/B(ep,k)) 1 



1 



-7k 



-1/2 



-1/2 



(27) 



This is the most general expression of finite temperature 
momentum distribution in the MI state in the large no 
limit. 

We find that above result can be also easily applied 
to the situation at critical point (N3 of Fig. [J), where 
u = 1 and ^ = U/2 in the large no limit. In the long 
wavelength limit, we have 7k 1 — ^a^k^, and therefore 



,N3 



coth 



U 



k a 



upto some overall density shit, Sn. When temperature 
is lowered, the amplitude becomes larger and hence the 
sharp peak is more pronounced. Note that in the low 
filling factor limit, tiq ^ 1, the exact cancellation in above 
derivation will not exist, and therefore the momentum 
distribution function will be still finite at k = 0, leading 
to a very sharp Lorentzian type distribution with a large 
correlation length. 



VIII. EFFECTIVE THEORY AND THE 
SCALING LAW OF THE MOMENTUM 
DISTRIBUTION NEAR THE CRITICAL POINT 

The quantum critical theory for the momentum distri- 
bution near the SF-MI quantum critical point at integer 
filling is belong to the univerality class of the d+1 di- 
mensional XY model [2^ (here d is spatial dimensional- 
ity). Before performing a full calculation, several remarks 



have to be mentioned: (i) It is known that the superfluid 
phase has only two independent fluctuations and there- 
fore we can use XY model or quantum rotor model with 
degree of freedom TV = 2 to investigate the quantum crit- 
ical phenomnea near the SF-MI transition, (ii) Since the 
anamolous dimension, 77, is known to be just 0.03 for 
D — 2 dimension and has a small logarithmic correction 
for 3D system [2^, we can expect that this problem can 
be very well approximated by the large TV result, which 
is equilvalent to the meanfield approximation in the limit 
of large filling fraction (n 3> 1). (iii) It is well-known that 
for 2D XY model, there is a classical Kotcrlitz-Thouless 
transition, which has a power-law decay of the correlation 
function in low temperature limit. The quantum critical 
phenomena of such quasi-long-ranged order is highly non- 
trivial and therefore will not be included in our following 
analysis. 

Based on above assumptions, we can consider the 
quantum field theory for large TV limit and use the value 
of saddle point to write down the meanfield order param- 
eters. Following the notation in Ref. [1^, the resulting 
fluctuations of action above such meanfield state can be 
written as the following XY model (or say the quantum 
rotor model of TV = 2 degree of freedom): 



TV ■«-^ f t^l'^ 
S = j dv dr {{drn^f ^ c^i^n^f 



~ 2 21 
-m nJ 



(29) 



where n^ for a = 1, 2 is the component of quantum fluc- 
tuations. (The original Hubbard model has U{1) sym- 
metry, so that a complex order parameter with two in- 
depedent components of quantum fluctuations are ex- 
pected.) Here c is the renormalized phonon velocity, g 
is the renormalized interaction, and m is the renormal- 
ized mass term, which changes sign at the SF-MI phase 
transition boundary [2^l23j . Note that the actual values 
of g, c and rh have to be calculated from a microscopic 
theory after integrating over the high energy modes and 
therefore can be different from the meanfleld results as we 
derived above within various regimes. However, since we 
are more interested in the universal form of momentum 
distribution rather than the quantitative values in this 
section, we will just keep them as a fitting parameters 
that might be able to be measured from the experimen- 
tal data. Besides, in this paper we are interested in the 
results in the low energy (or long wavelength) limit, and 
therefore we can neglect the lattice potential 

Here we first investigate the momentum distribution 
function in the normal state, which is a paramagnetic 
state in the spin (i.e. quantum rotor model with TV = 2) 
language. Using the notation in Ref. [1^, in the large TV 
theory, the correlation function between different com- 
ponents of the quantum rotor is given by the following 
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fomula: 



Ca/3(r,r) 

Xa0{k,U)n) 



(n„(r,r)n^(0,0)) (30) 

X{k,UJn)Sa,l3 (31) 

'^•^/^ (32) 



c^fc^ — (cj + i(5)^ + ' 



which is the single particle spectral function, and has 
been assumed to be isotropic in the low momentum 
regime (i.e. when the fluctuation length scale is much 
larger than lattice constant). As a result, the suscepti- 
bility X above depends only on the value of momentum, 
k = |k|. uJn — 2TrnT is the Matsubara frequencies, and 
we have let fc^ and h to be unit. In the last line of Eq. 
([5T|) . we have used the fact that the correlation function 
of a paramagnetic state (i.e. normal state) should have 
the diagonal part only in the spin component. 

The momentum distribution, n{k), is obtained from 
the equal time correlation and therefore 



Kk) = T^x{k,i^n) 



-cg/N 



c^k^ +10^,+ rh? 



(33) 



Note that different behaviors of n{k) in the different 
places of the phase diagram (Fig. [T](a)) just related 
to how the summation is calculated and how the mass 
term, m, changes in different regime near the critical 
point. Using the following identiy: Yln=-oo = 
(tt/x) coth(a;7r), we can evaluate the spectral function and 
obtain: 



n{k) 



cgl2N 



: coth 



\/(?k'^ + 
2T 



(34) 



When considering large T limit and finite mass term 
(to 7^ 0), we use coth(a;) l/x as a; — > 0, we find 



n(k) 



igT/N _ cgT/mN 

2 - TTeffc^' 



(35) 



which is the Lorenzian function one usually discussed 
in high temeprature regime, — c/fh is the coherent 
length. However, if we are looking at low temperature 
and/or large momentum regime, Eq. (|34p gives 



n(fc) 



~c~g/2N 



(36) 



as fc ^ to/c = ^jT^. This result implies the different mo- 
mentum dependence can be obtained when the momen- 
tum is larger or smaller than the inverse of correlation 
length, which is the length scale indicating the coherence 
of the entire system. If is large (i.e. lower tempera- 
ture and near the SF transition), most of the momentum 



distribution will be dominated by the quantum fluctua- 
tions, oc . On the other hand, if is small (i.e. away 
from the SF transition), the most momentum distribu- 
tion will be dominated by thermal fluctuations, which 
gives a Lorentzian type distribution. 

Finally, if we consider the momentum distribution 
along the the phase transition boundary or inside the 
ordered phase at finite temperature, i.e. m = 0, we can 
find that at finite T: 



n{k) 



~g/2N f£\k\ 
gT 



Nck^ 
9 

2N\k\ 



T:^dk 



T <^ck 



(37) 



This shows two different behavior for large T (or small 
k) and small T (or large k regime). This result agree 
with our statement from the three state model that ther- 
mal fluctuation is dominant and diverges as k~'^ in the 
small momentum regime, while the quantum fluctuation 
is dominant and scales as in the large momentum 

regime. 



IX. EXPERIMENTAL IMPLICATION 

The TOF image is known to be the most important 
measurement in the systems of ultracold atoms. The 
long time of flight image is usually believed to be equili- 
vant to the momentum distribution of the system inside 
the optical lattice. Our results show that how the dif- 
ferent scaling behavior and temperature/interaction de- 
pendence of such distribution of non-condensate particle 
can be distinguished and identified. This is of particu- 
lar importance when near the quantum critical point of 
the SF-MI transition. From experimental point of view, 
there are several methods to observe these different be- 
haviors in the current experimental setup. 

First of all, one can systematically investigate the 
TOF image for different initial temperature, and compare 
their TOF image both in the small and large momentum 
regime. According to our results, one should find that 
the TOF in the large momentum regime has little tem- 
perature dependence since it is mostly contributed by 
quantum fluctuations. On the other hand, after extract 
out the temperature dependent part of the TOF, one can 
use our analytical form of the thermal excitation (Eqs. 
PU)) and (|13p ) to obtain the temperature and condensate 
density inside the optical lattice. 

In recent experiments in Stanford's group [t'I, the TOF 
image of a condensate confined in the ID optical lattice 
has been fitted and investigated by using three Gaussian- 
type function: one for condensate, one for thermal atoms, 
and one for quantum depletion. Although the system 
setup and parameter regimes in their experiment are dif- 
ferent from the regimes we discuss in this paper, the clear 
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evidence for different temperature behaviors in the non- 
condensate particles still qualitatively supports our pre- 
diction: the quantum depleted particles dominates the 
large momentum distribution, while the thermal particles 
dominates the small momentum regime. Further inves- 
tigation on the TOF image can be very interesting and 
useful for the future development. 

Secondly, one can also separate the contribution of the 
thermal excitations from the quantum deplation by pro- 
viding a potential gradient just before releasing the opti- 
cal lattice. It is found [l^l that due to the different linear 
response of these two non-condensate particles, the cen- 
ter of thermal particle distribution can be away from the 
quantum deplated particles, making the direct measure- 
ment of these two distribution functions available. This 
approach opens the possibility to quantitatively study 
the condensate and non-condensate profile near the quan- 
tum critical point. 



X. SUMMARY 

In this paper, we have explicitly shown that the 
momentum distribution of the quantum deplated non- 
condensate particles can be very different from the ther- 
mal excited noncondensate particles, especially when 
near the quantum critical point of SF-MI phase tran- 
sition. Our results is also consistent with the general 
scaling properties derived from N = 2 quantum rotor 
model (i.e. XY model). Our analytical results can be 
used to provide a unique method to determine the tem- 
perature and condensate fraction inside the optical lat- 
tice. Futher extension to the study of noise correlation 
function and/or other physical properties can be also ex- 
pected in the same frame work. 
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(iii) WR(r) is the Wannier function, centered at a given 
lattice position, R. w-[i{r) is a Fourier transform of ^k(r) 
with respect to the lattice momentum, k. (iv) Finally, 
WR(q) is the Fourier transform of WR(r) with respect to 
the real space coordinate, r. Their relationship are given 
as following (L is the number of total lattice point along 
each direction and a is the lattice constant, ' means 
a summation over the first BZ): 



*k(q) 
*k(r) 

«'k(r) 

WR{r) 

WR(q) 
^'k(q) 



J dre-*i ""«'k(r) 

k 

I (0^^^''-^^-(^) 

k 



(Al) 
(A2) 

(A3) 
(A4) 

(A5) 
(A6) 

(A7) 
(A8) 



Note that the Fourier and the inverse Fourier transform 
convention is different for real space momentum and lat- 
tice momentum. The amplitude of a wavefunction can 
be shown to be normalized: 



1 = / dr|*k(r)|^ = 



l*k(q) 



(2^)3 ' 

drMv^^ I ^|^;R(q)|2 (A9) 



APPENDIX A: MOMENTUM DISTRIBUTION IN 
AN OPTICAL LATTICE AND IN FREE SPACE 

1. Definition of wavefunction and tiieir Fourier 
transform 

We first identify the following four kinds of single par- 
ticle wavefunctions in the optical lattice and clarify their 
relationship with each other, (i) vE'k(r) is the Bloch wave- 
function in real space at a lattice quantum number, k, 
in the first Brillouin Zone (BZ). (ii) 5'k(q) is its Fourier 
transform with respect to the real space coordinate, r. 



2. Momentum distribution 

Now we need to identify the field operators based on 
the following three different basis: (i) lattice site ba- 
sis (R) with Wannier function (wR,{r)) as an eigenstate, 
(ii) Bloch momentum basis (k) with Bloch wavefunction 
(^k(r)) as an eigenstate, and (3) free space momentum 
basis (q) with plane waves (^^j and V is the volume of 
the whole system of free space) as an eigenstate. Their 
fields operators are defined to be cr, ak and 6q, respec- 
tively, and they are related to the single particle field 
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operator, ^'(r), to be 
#(r) 



y^aRWR(r) 

R 



V 



(AlO) 



Therefore, these field operators can be transformed with 
each other as following: 



1 



Ok 



L3/2 ^ 
R 



Org 



(All) 
(A12) 



As a result, the momentum distribution in real space 
{N{q)) can be connected to the momentum distribution 
in lattice {n(k)) by 



N(q) ^ (6t 6 ) 



where we have used (ajj-ak' 



q q' 

k 



(A13) 



7i(k)Sk,k'. 

It is known that when the lattice strength is strong, 
we can approximate the Wannier function of the lowest 
subband by a Gaussian-type function: 



WR(r) 



1 



,-i|r-R|V-^ 



7^3/4^3/2 

where a is the Gaussian width. Therefore, we have 

2 



(A14) 



R 



where 



/(q) 



|w(q) 



: /(k-q)Kq)p 

(27r)3 

-3/2' 



(A15) 



E e'"" ^ 

R 



n 



sin(L(7cta/2) 



sin((7ja/2) 



(A16) 



(A17) 



In the limit of large system size {L — > oo), /(q) becomes 
a periodic function of delta-fuction peaks at each recip- 
rocal lattice, Gn = (27rn/a), where n = (ni,n2,n3) is a 
vector of integer numbers. Within such approximation, 
the momentum distribution in free space can have a sim- 
ple relationship with the momentum distribution in the 
optical lattice: 

7V(q) « EE'N'')Hq)l''^(k-q-G„)](A18) 



However, such simplified result is not correct if the sys- 
tem of finite size, where /(q) is not a delta function, but 
just a sharp function at Gn. 



3. Momentum distribution after column 
integration 

The TOF image taking after a long expansion time 
can be understood as the momentum distribution after 
column integration along a certain direction (say z axis) . 
Therefore, it is instructive to study the momentum dis- 
tribution in free space after column integration along the 
z axis. From a general result of Eqs. (jA15P - (jA17|) . we 
have 



N 



(qx) = !:'»«/ If |*I.(<1) 



2tt 



E 

R,R' 

lE'"(k)|^k,(q^) 



L3 ^3/2 

^»(k-q)-(R-R')-q''^' 



V 

k 

where (q_L ) is the Fourier transform of the 2D Bloch 
function in 2D lattice, similar to the 3D case defined in 
Eqs. dMI, (Em, (TCTel) . and (ElTl) . 

To evaluate the summation over and R'^, we note 
that the first order cotribution comes from the case Rz = 
while the second order contribution comes from \ Rz — 
R'zl = a, which is exponentially smaller the first order as 
a S> <T near the SF-MI transition regime. As a result, it 
will be a good approximation to keep the leading order 
contribution only and hence the momentum distribution 
in free space after column integration becomes 

A^j_(q_L) oc E'"-L(k_L) *k^(q_L) 

k_L 

cx ni(ki)|w;(qi)p, (A20) 

if we are interested in the long wavelength limit (q ^ 0) 
and approximate the function /(qj_) oc S{q± — k_L) in the 
limit of infinite number of lattice (L — > oo). Here 



tE'^w 



(A21) 



is the momentum distribution in optical lattice after col- 
umn integration. 
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APPENDIX B: FINITE TEMPERATURE 
MOMENTUM DISTRIBUTION IN SUPERFLUID 
PHASE 

1. Canonical transformation 

In order to clarify the complicated process, here we list 
all the transformation of quasi-particles used in the SF 
phase within the three-state model. Within the three- 
state model, Altman et al. [l6j introduced a new oper- 

ator, t\^\Q) = (aRj |0) = |no -I- a)^ to 

create a state of no + « particles (a = ±1 and 0) at site 
R. As a result, one can show that within this trucated 
space, any local operator, say the bosonic field opera- 
tor Or, and the particle number operator, a^Op^, can be 
respectively expressed to be 

an = Vno + l4,R^i,R + \/^^-i,R^o,R (Bl) 

ak^R = ('^0 + l)4,R^l,R + '^O^o.R^O^R 

+ {no - l)tli3i_i.R. (B2) 



We note that in the three-state presentation, the number 
operator, a^^Oj,^, is not a simple product of the result 
for Or and a^. This can be easily understood from the 
fact that when applying a^^aj^ on a state \no — 1)r, the 
particle number can be virtually reduced to no ~ 2 by a-^ 
first then becomes its original value by a^. such virtual 
process is not included in the three-state representation 
of flR. Such trivial notational problem occurs only in the 
situation when two or more field operators applying on 
the same lattice site. 



Within the meanfield approximation, we can assume 
that the SF ground state is a coherent state com- 
posed by different number of states at each site, i.e. 
I*mf) = IIr&UIO)' ^liere bl-^ = cos(0/2)4_j^ + 
sm{9/ 2) {cos xtl r + sinx^Li r) where 9 and x are two 
variational parameters [IB] . When considering the quasi- 
particle excitation, Altman et al. further generalize 
above transformation to the following: 



''O.R 

if 



-1,R 



cos(6'/2) 



■sin(6l/2) 







i(0/2)cosx cos(0/2)cosx — sinx 
sin(0/2)sinx cos(0/2)sinx cosx 



Uqo Uqi Uq2 

Uio Un Ui2 
U20 U21 U22 



4,R 



= U-b 



"0,R 

^Ir 



^2,R 



(B3) 



where U is the transform matrix, and 6r — 
[^0 R' ^1 R' ^2 r]"^ vector of operators with 61/2, r 

being the field operators for the excitations above the 
condensate state, foj ^. 

When in the SF state, one can approximate the sin- 



gle particle ground state as a c-number, and apply the 
generalized bogoliubov tranformation to diagonalize the 
system Hamiltonian upto the qudratic order of excita- 
tions. Such generalized Bogoliubov transformation can 
be written to be [IB] 



t 

k 



iVn(k) iVi2(k) Pn(k) PM 

Ni2{k) iV22(k) Pi2(k) F22(k) 

Pii(-k) Fi2(-k) iVn(-k) iVi2(-k) 

Pi2(-k) F22(-k) iVi2(-k) iV22(-k) 



(B4) 



where the definition of the excited states, 6k and /3k are we will use numerical method to evaluate them directly, 
clearly shown as above. Al(k) is the 4x4 transform 
matrix, which in principle can be derived analytically. 
However, since the full analytic expression is too compli- 
cated to be understood in a simple physical picture, here 
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2. Momentum distribution of noncondensate We first transform it into the ^-operators: 

particles 

Our goal is to calculate the momentum distribution in 
optical lattice: (aj^Ok), by usng above transformations. 



R,R' 

+ y/no{no + 1) (^(4,R^o,R^-i,R'*o,R') + (4,r*-i,r4,r'*i,R') , 
The first term, defined to be Ti upto the prefactor, tiq + 1, gives 

1 ^ e''^ (^-^')(4_i,io.R4,R'ii.R') =^'<5k,oC/i'oC/oo 



(B5) 



R,R' 



R,R' 



R 

ik-(R-R') 



o 

;:^l'o((C^Ol6l,R + C/02623)(f^01^lR' + ^^024^')) + U^oiiUllbl^ + C/i24,r)(C^11^1,R' + t^l2&2,R')> 

E e^'"^''"'''^ [{{Uoib,,n + C/o262,r)(C/ii6i,r, + C/i262,r')) + ((C^n^^U + C^i24,r)(C^oi6U- + C/o24.r')) 



R.R' 



(B6) 



where we have used the approximation: feo.R ^ 1 — 
l4 R^i R ~ i4 R^2 R eliminate 69. r. to the quadratic 
order of small fluctuation, 61/2, r- The first two terms are 
proportional to the system size, L^, and are nonzero only 
at k = 0, indicating a contribution from the condensate 
particles. The last two terms, however, are the noncon- 
densate part. Since in our paper, we are interested in 
the momentum distribution at finite lattice momentum 



J 



k, we will concentrate on the later part for the numerical 
calculation of momentum distribution. 

Below we define the momentum distribution of the 
noncondensate particle from the ith term of Eq. (|B5|) 
to be Ti {i = 1,2,3,4) (not including the prefactor, uq, 
no+1, or -y/ no{nQ + 1)), and these terms can be expressed 
in the momentum space to be 



Ti = 



To = 



T. = 



bU 
bL 



2,-k 



1. k 

2, k 



^l.-k 
5'2.-k 



"^00 "^11 

UiqUoqUqiUii 
UiqUqoUq2Uii 

'^20'-'01 
t^20^02C^01 

U2oUoaU2iUoi 
U2QU00U22U01 

U0QU20UIIU01 
U00U20U12U01 

UioU2oU§i 
UioU2aUoiUo2 



U§oUiiUi2 UioUooUiiUoi 

^00^12 UiqUoqUi2Uoi 
U10UQ0U02U12 UiqUoiUo2 



U20UQ2U01 

"^20 '-'02 
U20U0QU21U02 
U20U00U22U02 



U20U00U01U21 
U20U00UQ2U21 



U^nU2lU22 



U00U20U11U02 U^oUllU21 

UQoU2aUi2Uo2 UqqUi2U21 

UioU2oUa2Uoi UioUoaUQiU2i 

UiaU2oUQ2 UioUoaUQ2U2i 



U10U00U11U02 
U10UQ0U12U02 

UiqUoiUo2 
'^10'^02 

U20U00U01U22 
U20U0QUQ2U22 
UqqU2iU22 
'^00'-'22 

U§oUllU22 

U§qU 12U22 
UioUoQUaiU22 
U10U0QU02U22 





' b,M - 








bU 




UuJ 








&2.k 




bU 




L4,-kJ 








&2.k 




bl.-k 




L4,-kJ 



(B7) 

(B8) 

(B9) 
(BIO) 
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and the forth term (T4) is the exactly the complex con- 
jugate of the thnd term (T3). Therefore after mclud- 
ing the proper prefactor in Eq. (jBSp . the total effect of 
non-condensate particle can be obtained by adding them 
together to be 

(aj^ak)k5<io = {m + l)Ti(k) -|- noT2{k) 

+ V"oK + l)(T3(k) + T4(k)) 

= bl ■ r(k) • 6k, (Bii) 

where 6k = [&i k' ^2.k' ^1 -k' ^2 -k]^ ^^'^ '^0^) is the re- 
sult of the final 4x4 matrix {not the same as the 3x3 
transform matrix, U, defined in Eq. (IB3P ). 

In the final step, we have to transform the operator of 
^1/2, ±k to /3s/m,±k (as shown in Eq. (jB4p ) to calculate 
the expectation value. The obtain momentum distribu- 
tion from the noncondensate particles can therefore be 
written to be 

{aia^,)k^a = bl ■ T{k) ■ 6k 

- /3^X(k)t.T(k)•M(k)Ac 

= Pi ■ S{k) ■ /3k, (B12) 

where the new matrix 5(k) = A^(k)^^-T(k)-A^(k). Using 
the fact that {Pl/^^^Ps/mM) = fsi^s/m,)^) for the Bose- 



Einstein distribution of particle/hole excitation. There- 
fore, the final result of the noncondensate particle mo- 
mentum distribution can be expressed by the matrix el- 
ements of matrix 5(k); 

(4«k>k#0 = (5ll(k)2 + 533(k)2)/(e,,k) 

+ (522(k)2 -I- 544(k)2)/(e™,k) 

-K533(k)' + 544(k)^ (B13) 

where the last two terms are the contribution of quantum 
deplated particles. 



APPENDIX C: FINITE TEMPERATURE 
MOMENTUM DISTRIBUTION IN MOTT 
INSULATOR PHASE 

1. General expression in three-state model 

In the Mott insulator phase, the meanfield ground 
state has one particle persite. Therefore it is important 
to separate the contribution of onsite density operator, 
^R^R' from other non-local density operators. In order 
words, when calculating the momentum distribution, we 
have 



Fl.i,Fl.2 R. R.i^R.2 

E^ [("0 + l)4,R^l,R + '^o4,rVr + ('^O - l)<ll,Ri-l,R] ) 
R 

^Jl E^f^'^0 + l4,R^0,R + \/^4,R^-l.R. V'^O + l4,R*l,R + \/^i-l,R*O.R ) 
R 

+ ^3 E ^[^^^O + l4,Ri^0,Ri + \/^4,Ri^-l,Ri n/"0 + l4.R2*l,R.2 + V^*-132^0,R2 



R-i ,R,2 



„ik-(Ri--R2) 



(CI) 



where in the last line, for the convenience of later Fourier 
tranform into momentum space, we have added the term 
of Ri = R2 and substracted it in the second term. 

Before transforming to the momentum space represen- 
tation for the calculation of quasi-particle distribution, 
we can simplify above results further by using the num- 



ber constrain, ^L-i 4o+a,R^«o+a,R ^ ^ ^o replace the 
field operator, io,R and to the quadratic order of 
small fluctuation, After some straightforward calculation, 
it is easy to show that 



aj^Ok) 



no 



Jji E + l)(4,kii,k> + "o(i-i,ki-i,k> + v/"o(«o + 1) ((i-i.-k^i.k) + (4,ki-i,-k>) 



+ ("0 + l)(4,k^i,k> + "o(^-i,ki-i.k> + v/"o(no + 1) f(i-i,-k<i,k) + (4,k^-i,-k> 



(C2) 
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where in the last hne we apply Fourier transofrm, ta.k = 
jjj2j2R^a,Re~^^'^, to the momentum space represen- 
tation. Note that the first line of the last equation is 
independent of momentum k, and therefore can be un- 
derstood as a constant shift of the momentum distribu- 
tion. 

I 

(4,kii,k> 

where fsiep/hM) = (/^p//, -k/'p/Zi -k) is the Bose-Einstein 
distribution of particle/hole excitations. 

Therefore we can obtain the final expression of the mo- 
mentum distribution, n(k) = {al^ay^), to be 

n(k) = ^Sn + fBiep,k){V^^^A{k)-^Bik)f 

+ (1 + fBieh.k)) (V^^^Bik) - V^A(k))' , 

(C6) 

where the density shift, dn is defined by 

k 

+ {VmTTB{k) - V^A(k)) (1 + fBieh.k))] - "0 

(C7) 

to ensure the conservation of total number of particles. 



2. Analytical calculation of chemical potential 

In order to evalulate the momentum distribution in the 
normal state regime, the first thing is to determine the 
chemical potential by fixing the total number of parti- 
cles. We can use the general expression of momentum 
distribution, i.e. (aRaR.) in Eq. (jCl[) . and obtain the 
total number of particles to be 

N = E'(4«k>=E(«R«R) 

k R 

= L^np + y^(4.R^i,R - ^-i.r^-i,r) 

R 

k 

Now if we tranfer the i— operator to the quasi-particle 
operators, /3p/h,k, in the MI state via Eq. ([5]), the con- 
servation of total number of particle can be expressed to 



We now calculate the expectation value for each term 
by transforming to the quasi-particle, /3k, field via Eq. 
([5]) and neglecting those particle non-conserving terms 

(for example, (/3p,k/3p,k) = (/3l.k/?I,k> = ^ 0)- 

We have 



(C3) 
(C4) 
(C5) 

I 

be 

= E W-^('^)')[/s(^P^k)-/BKk)] 
k 

= J2lfBi^pM)-fBiehM)], (C9) 

k 

which is nothing but the difference of particle excita- 
tions and hole excitations. Below we will investigate an- 
alytically the chemical potential in different parameter 
regimes of interest. 



a. In the low temperature limit, U ^T,J 

We first consider the situation when U ^ T,J, so that 
we can use large U expansion for the particle and hole 
excitation spectrum and obtain 

ep,k = ^ - (3^+ {no + l)eo(k)) 

- "°^"°+J^^°(^)' +0((J/L/f) (CIO) 
e;i,k = ^ - {-5fi + nQeo{k)) 

-!^^^I^^l±^}^ + o{iJ/un (cii) 

Since U :s> T, we can expand the Bose-Eistein distribu- 
tion fuction as 

1 

~ g(;7/2-5e(k)T5M)/T „ I 

^ p-(c//2-5.(k)T5p)/T _^ oie-^/^),{G12) 

where ^e(k) = (hq + l)eo(k) for particle excitation and 
6e{k) = noeo(k) for hole excitation. Now using the fact 
that /^^ ge^™s(fc) ^ j^^^)^ ^Yiere Io{x) is the modified 



- A{k)^(3l_^P^^_^)+B{k)'{f3,,^^Pl^)^A{kffB{ep^^^ 
= B(k)2/B(ep.k) + A(k)2(l + /s(e,.,k)), 

= (^-i,-k4,k>* = -^(k)i?(k)/s(ep,k) - A(k)i?(k)(l -H/s(e„,k)), 



17 



Bessel function of the first kind, we can integrate out the 
momentum integral of Eq. (jC9p and obtain 



-y 

k 



and 
1 



7i3(ep(k)) 



g-t//2Tg5M/T 



Iq 



2{no + 1)J 
T 



1 3 



(C13) 



(C14) 



Here we have neglected the higher order terms, which is 
proportional to e~^/"^ <C 1. As a result, the chemical 
potential can be calculated to be (from Eq. jCU) 




(C15) 



which is good for all the temperature range if only U 3> 
J, T. If consider intermediate temperature regime, U 3> 
T ^ J, which is the most relevant regime in the present 
experiment, the leading order result gives 



J 4^ 



(C16) 



where z = 6 in 3D and z = 4 for 2D systems. 

On the other hand, in the limit of small T (i.e. T J), 
we use Io{x) — \/^e^ and obtain 



no 



no 



3/4 



or 



5/i 



zJ 



8T fno + 

— log 

z \ no 



This is consistent with the results at T = 0. 



(C17) 



(C18) 



b. In high temperature interaction limit, T ^ U ^ J 

In the limit of high temperature, we can expand the 
excitation modes and Bose-einstein distribution to the 
leading order terms of small J/T and J/U. we have 



(5e(k) ± 6)1 



+ (eC//2T_i)2 (^gC//2T_i )\ T ) ^^^^^ 

Therefore within the same order, the total number of particle and hole excitations shown in Eq. (jC9p can be calculated 
to be (note that eo(k) = 2J^^ cos(fcQa) for 3D optical lattice.) 



^ 2^ JbK^pW) - ^ui2T _ + (gC//2T _ 1)2 y + {e^/'^T _ ly y(,U/2T _i 2^ /Vr2+ ^•'"0 + J-2 j 

(C20) 

^ 2^ /B(e/i(k)) - g[//2T _ 2 ^ (e^/2^ - 1)2 T (e^^/^^ - 1)^ Ve^/^?' - 1 2*^ j y2 + ^"o ^2 j -(.^^^^ 

I 



The chemical potential then can be determined by the 
following colse form to the second order of J/T 

in the high temperature and strong interaction limit, U ~ 
T> J. 



3. Analytical calculation of momentum distribution 

a. In the low temperature limit, U ^ T,J 

In order to study the momentum distribution shown in 
Eq. (jC6[) . we can start from the following leading order 
expansion in the strong interaction limit {U ^ T, J) by 
us ing Eq . ©: (V^o + lS_(k) - ^^^^^^(k))^ ~ no, and 
(V"-o + l^(k) - ^/n^BCk))'^ ^ 1 + no. Now, after in- 
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eluding the large U expansion of particle-hole excitation The momentum distribution then becomes 
shown in Eq. (|C12p . the density shift can be expressed 
to the leading order to be 



6n 



-1 ' + l)e-^/2Tg(„o + l)eo(k)/TgWT 

k 

+no (l + e-t^/2^e"°^°('')/^e-*^/^)] - no 



(C23) 



i(k) = {V^^^A{k)-^B{-k)ffB{ep,^) + {V^^^B{-k)-^A{-k))^{l + fB{ehM))-Sn 



no + 2no(l + no)^^ + (1 + no) 



2(710 + 1)./ ^ 

T ^ 



g~C//2Tg(no + l)eo(k)/T 



,(C24) 



to the leading order of large U, i.e. the next higher 
order is proportional to {J /U) e~^/'^'^ . It is easy to 
show that the total number of particles is conserved, 
i.e. "/T'X^k'^C*) ^ ""0, and therefore above expression 
is a close form for the momentum distribution in large 
U limit. It is also straightforward to calculate the higher 
order terms but we will not show them here. From the 
momentum distribution in the optical lattice, we can also 
integrate out the z-component of the lattice momentum 
and obtain 



where eo(k_L) — 2J(cos(kxa) + cos(A;j,a)) and 



= no — ^n + 2no(l + no) 



eo(k±) 
U 



Gi = 



Go 



T ( 2naJ \S/2 
t//2T ^0l T (C26) 

/o(^)i/2 



(C25) 



6. temperature limit, T ^ U ^ J 



For the high temperature limit, we need to expand the 
whole formular in terms of J /T and J /U at the same 
time. To the leading order, we thus have 



i(k) = (V^^;;TTA(k)-y^B(k))Vs(ep,k) 

eo(k)\ 1 



1 + no + 2?io(l + no) 
no + 2no(l + no) 



U 

£o(k) 

U 



eC//2T _ 1 



1 - 



1 



et//2T _ 1 



no + 27io(l + no) 



eo(k) 4no(l +no) eo(k) 



U 



JJ/2T _ I 



u 
r 



lB(k) - V^^A(k))2(l + jBiehM)) - 5n 



pU/2T 

1 



T 

noeo(k) - 5^ 



et//2T _ I T 
(1 + 2no + 2n2)e^/2^ eo(k) 



5n 



[eU/2T _ ly 



T 



(C27) 



where the correction of chemical potential is the second last line above, we also used the fact that 5n is deter- 
order effect (Eq. (|C22|) ') and has been neglected. In the 
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mined by requiring X]k''^(^) = "-o *° leading or- 
der, and therefore obtained 5n = -grf^r^- Therefore, 
after integration over the component, we have 



where 



Hi = 



n_L(k_L) = riQ + i?ieo(k_L)/J, 
2rto(l + no) J 



U 



1 



pC//2T _ I 



(1 + 2no + 2n2)e^/2T j 



(eC//2T „ 1)2 



T' 



(C28) 



(C29) 
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